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Abstract. In a recent paper, Clusel and Fortin [J. Phys. A.: Math. Gen. 39 
(2006) 995] presented an analytical study of a first-order transition induced by an 
inhomogeneous boundary magnetic field in the two-dimensional Ising model. They 
identified the transition that separates the regime where the interface is localized 
near the boundary from the one where it is propagating inside the bulk. Inspired by 
these results, we measured the interface tension by using multimagnetic simulations 
combined with parallel tempering to determine the phase transition and the location of 
the interface. Our results are in very good agreement with the theoretical predictions. 
Furthermore, we studied the spin-spin correlation function for which no analytical 
results are available. 
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1. Introduction 

Wetting transitions arc phase transitions in the surface layer of bulk systems which are 
induced by symmetry-breaking surface fields [1, 2]. The Ising model with a boundary 
magnetic field is a simple model for such a wetting problem, because Ising ferromagnets 
have the same critical behaviour as the analogous case of gas-fluid transitions, as has 
been pointed out by Nakanishi and Fisher [3]. The use of the Ising model with short 
range interactions for wetting studies has not only the advantage that one can use 
all the advanced simulation techniques which have been developed in the past years. 
Especially in two dimensions (2D), there are also a lot of theoretical results available 
for comparison. 

The Ising model with a uniform boundary magnetic field on one side of a square 
lattice has been completely solved by McCoy and Wu [4], whereas the Ising model with 
a uniform bulk field can only be solved at the critical temperature [5]. For situations 
with fixed boundary spins or equivalently infinite boundary magnetic fields [6] , or finite 
boundary magnetic fields [7] some exact results have also been found. In a recent 
paper, Cluscl and Fortin [8] presented an alternative method to that developed by 
McCoy and Wu for obtaining some exact results for the 2D Ising model with a general 
boundary magnetic field and for finite-size systems. Their method is based on the 
fermion representation of the Ising model using a Grassmann algebra. They applied 
this method to study the first-order transition induced by an inhomogeneous boundary 
magnetic field in the 2D Ising model [9]. To be more precise, the boundary magnetic 
field acts on the x — 1 column of spins, being positive in the lower and negative in the 
upper halve. By taking the thermodynamic limit exactly for a given geometry of the 
lattice, they obtained a simple equation for the transition line and also a threshold for 
the aspect ratio = L,j./Ly = 1/4, where this line moves into the complex plane. This 
vanishing of the transition line indicates the crossover from ID behaviour for <^ Ly 
to 2D behaviour at large which is reflected in the behaviour of the boundary spin-spin 
correlation function. 

The aim of this work is to check some of the predictions by carrying out Monte 
Carlo simulations of this model and to extend the results to parameter ranges and 
for observables where analytic solutions cannot be obtained. The rest of the paper 
is organized as follows. In Section 2 we give the definition of the model and briefiy 
summarize the theoretical predictions. A description of the employed simulation 
techniques and the results of our Monte Carlo simulations are presented in Section 3, 
and concluding remarks can be found in Section 4. 
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Figure 1. The phase diagram at zero temperature as a function of the aspect ratio 
Q = Lx/Ly and the boundary magnetic field H. 



2. Model and Theoretical Predictions 

We consider a 2D Ising model with a non- homogeneous magnetic field hy located on one 
boundary of the system. The Hamiltonian is given by 

'H — — J ^ {cTxyCTx+ly + CFxy(^xy+l) — ^V^^V ' (-'-) 
x,y=l y=l 

with free boundaries in the x-direction and periodic boundary conditions in the y- 
direction. To compare our results with the theoretical predictions of Clusel and 
Fortin [9], we consider the same profile of the boundary magnetic field acting on 
the X = 1 column of spins: hy = H for y = l,...,Ly/2 and hy = —H for 
y = Ly/2 + l,...,Ly, with H >0. 

In the limit of zero temperature, by using simple energetic arguments, Clusel and 
Fortin [9] showed that for small H all spins are ahgned in one direction as in the bulk 
case for H — 0, see also Figure 1. With increasing H, however, depending on the aspect 
ratio ( = Lx/Ly two different interfaces can be formed. If ^ > = 1/4 the interface 
is localized near the boundary, whereas for C < Cs the interface is propagating inside 
the bulk. The critical ratio (s marks the crossover from a ID behaviour for Lx ^ Ly 
towards a 2D behaviour at large (. For ( < (^s and non-zero temperatures T > 0, 
with the abbreviations t = tanh( J/Zc^T) and u = ta.nh.{H /ksT), the equation for the 
first-order transition line in the {t, M)-plane turns out to be a quadratic equation in u'^ 
[9]: 

2t(l + ^;(4C))«' + (1 + t^) (l - 2tv{4C) - t^)u' + 2(^;(4C) - l)t^ = 0, (2) 



^(40 = cosh 



, U(i + *), 
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Figure 2. The phase diagram for a system with C = 0.2. The thick hne shows 
the first-order transition given by Clusel and Fortin [9] and the thin vertical hne 
indicates the second-order bulk phase transition. The thin dashed line indicates a 
simple approximation of the first-order line, H w 2^^^^, where cr{T) and TOo(T) are 
the interface tension and the spontaneous magnetization of the pure 2D Ising model, 
respectively. The sketches of the spin configurations illustrade the interface location 
in the three phases and the double-headed arrows show the parameters of the Monte 
Carlo simulations. 

In Figure 2 we show the phase diagram for a system with aspect ratio ( = 0.2 (and 
J = ks = 1). In the low-temperature regime we can approximate the above expression 
by comparing the energy of the interface with the energy induced by the magnetic field. 
This leads to if ~ 2a{T)(/mQ{T), where cr{T) and mo(T) are the known interface 
tension and the spontaneous magnetization of the pure 2D Ising model, respectively. 
This approximation reproduces the exact low-T expansion, H = AC, — 4(^Te~^/^, and 
works very well for T < 1 as one can see in Figure 2 (thin dashed line). Since o"(T) 
vanishes much faster than mo(T) as T — Tc, also this point is reproduced exactly, but 
the slope of the approximate transition line at does not diverge as for the exact 
solution. 

Due to the first-order transition induced by the inhomogeneous boundary magnetic 
field, the second-order phase transitions across the vertical line at T = Tc are transitions 
from a region where an interface in the bulk separates two ordered domains of opposite 
magnetization from a disordered regime above the transition temperature. Therefore, 
the system undergoes a transition without a change in the magnetization (m) which is 
zero in both phases, cf. Figure 2, but the width of the magnetization distribution does 
change. 
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3. Numerical Results 

Since we are primarily interested in the location of the interface induced by the 
boundary field, we first performed simulations at low temperatures to generate a well- 
defined interface. To overcome the slow dynamics at low temperatures we developed a 
combination of the multimagnetic algorithm with the parallel tempering method [10] 
for which we used two different schemes: In the first scheme, we kept the magnetic field 
value H fixed and simulated n = 32 replica of the system at different temperatures Tj. 
In the second scheme, we kept the temperature T fixed and used n = 32 different values 
of the magnetic field Hi. 

To construct the weight function for the multimagnetic part of the algorithm, we 
employed an accumulative recursion, described in detail in Refs. [10] and [11]. Statistical 
averages were taken over runs of 1 x 10^ Monte Carlo (MC) steps, where one MC step 
consists of one full multimagnetical lattice sweeps for all 32 replica and one attempted 
parallel tempering exchange of all adjacent replica. With this method we were able to 
study systems with = x L^^ = 50 to 5000 spins for aspect ratios C, — L^/Ly — 0.2, 
0.25 = C,s and 0.5, for further details sec Table 1. 

Let us first discuss the data obtained for the case = 0.2 < (s- For this value of 
the aspect ratio, the phase diagram as predicted by Clusel and Fortin [9] is shown in 
Figure 2. The thick line indicates the first-order transitions from the fully magnetized 
state with (m) > to the mixed state with an interface extending across the bulk. 
To check the nature of these transitions we measured the probability density of the 
magnetization at four points along the transition line. In the first two cases, we kept 
the boundary magnetic field constant {H — 0.4 and 0.5) and varied the temperature 



Table 1. Summary of simulation parameters (PT: parallel-tempering algorithm, SC: 
single-cluster update). 



c 


H 


T 


Lx 


X Ly 


method 


measurements 


0.2 


0.4 


1.6 - 2.2 


80 - 


- 2000 


PT 


1 X 10^ 


0.2 


0.5 


1.4 - 1.9 


80 - 


- 2000 


PT 


1 X 10^ 


0.2 


0.5 


2.1 - 2.3 


80 - 


180500 


SC 


1 X 10^ - 5 X 10^ 


0.2 


0.7-0.9 


1.0 


80 


- 640 


PT 


1 X 10^ 


0.2 


0.285 - 0.316 


2.0 


80 - 


- 2000 


PT 


1 X 10^ 


0.25 


0.4 


1.9-2.1 


256 


- 2500 


PT 


1 X 10^ 


0.25 


0.5 


1.5 - 1.9 


64 - 


- 2000 


PT 


1 X 10^ 


0.25 


0.5 


2.0 - 2.3 


64 - 


- 6400 


SC 


1 X 10^ - 5 X 10^ 


0.25 


0.7-0.8 


1.5 


64 - 


- 1600 


PT 


1 X 10^ 


0.5 


0.5 


2.2 - 2.35 


50 


- 3200 


PT 


1 X 10^ 


0.5 


0.9-1.1 


1.0 


50 - 


- 5000 


PT 


1 X 10^ 


0.5 


1.5 - 2.0 


1.0 


50 - 


- 3872 


PT 


1 X 10^ 
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Figure 3. Left plot: The probability density for the magnetization as a function of 
temperature at constant boundary magnetic field H = 0.5 for iV = 20 x 100. The 
vertical lines indicate the first-order transition temperature Tq w 1.5950 according to 
Equation (2) and the critical temperature Tc ~ 2.2692 of the bulk phase transition. 
Right plot: Histograms of the magnetization in the vicinity of To for various lattice 
sizes ranging from Ly = 20 to Lj, = 100. Here the temperatures Tq(L) are determined 
such that the peaks at m « ±too and m — are of equal height. 



to locate the transition point, and in the other two cases, we fixed the temperature 
(T = 1.0 and 2.0) and varied the boundary magnetic field. These points are indicated 
by the double-headed arrows in Figure 2. 

In the following we illustrate our procedure for obtaining the first-order transition 
point and the associated interface tension for the case of fixed H = 0.5. A level plot of 
the magnetization density m = (l/N) I]a;,j/'=i (^xy as a function of temperature is shown 
in Figure 3 (left). For each lattice size, a pseudo-transition point can be defined by 
varying the temperature until the peaks at m ~ ±77io and m = are of equal height, 
which can be achieved by histogram reweighting. The interface tension can then be 
estimated from [12] 



n = W7^^h^], (3) 




where P™^^ is the value of the peaks and P™™ denotes the minimum in between, see 
Figure 3 (right). The length of the interface is denoted by L, which is L = L^; in the 
case of C < Cs- 

The thus defined pseudo-transition temperatures To(L) approach the infinite- 
volume transition temperature Tq as l/L"^, and for the final estimate of = lim^^oo F[, 
we performed a fit according to 

b\n(L) 

P£ = p^ + - + _^. (4) 

At fixed T one proceeds analogously by varying the magnetic field H, i.e., the roles of T 
and H are just interchanged. For all four cuts at constant surface field or temperature we 
find a good agreement with the infinite- volume transition points derived from Equation 
(2) and a clearly nonzero interface tension, see Table 2. 
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Table 2. Results of the simulations close to the first-order transition line. For a given 
aspect ratio C we kept either the boundary magnetic field H or the temperature T 
fixed. The third column shows the measured transition points and the fourth column 
the exact infinite- volume values given by Equation (2). The fifth column contains our 
numerical estimates for the extrapolated interface tensions. 



c 


H 


To 


To (exact) 


Fs 


0.2 


0.4 


1.84(1) 


1.82252. . . 


0.18(1) 


0.2 


0.5 


1.60(1) 


1.59497... 


0.32(2) 


0.25 


0.5 


1.91(1) 


1.95845... 


0.09(1) 


C 


T 




Hq (exact) 


Fs 


0.2 


1.0 


0.72(1) 


0.702352. . . 


0.82(2) 


0.2 


2.0 


0.305(3) 


0.305928. . . 


0.12(1) 


0.25 


1.5 


0.73(1) 


0.762807. . . 


0.24(1) 



We also checked the critical behaviour along the line of second-order transitions 
at T = Tc = 2/log(l + \/2) K 2.2692. To this end we run at = 0.5 single-cluster 
simulations (suitably adapted to the surface field) for systems with iV = L^. x Lj^ = 4 x 20 
to 190 X 950 spins and performed a finite-size scaling (FSS) analysis to determine the 
transition point and some critical exponents. This particular value of the magnetic 
field has been chosen because of the relatively large temperature gap between the 
boundary field induced first-order transition and the bulk phase transition. Between 
each measurement we performed one sweep, which here consists of n single-cluster 
updates with n chosen such that n{\S\) N, where {\S\) is the average cluster size. 
For every run we generated 10^ sweeps, and recorded the time series of the energy 
density e = E/N and the magnetization density. Using these time series, we can 
compute the specific heat, C = N{{e'^) — (e)^)/T^, the (finite lattice) susceptibility, 
X — N{{m'^) — (m)^), and the Binder cumulant U — 1 — (m^)/3(m^)^ in the vicinity of 
the simulation point by reweighting. 

In this way we can use the maxima of the (finite lattice) susceptibility to detect 
the pseudo-critical points and can obtain an estimate for T^, from a linear least-square 
fit of their scaling behaviour, Tmax — Tg oc L~^^'^ = L^^, assuming thus the exact value 
u — 1 according to the universality class of the 2D Ising model. This leads to an 
estimate for the critical temperature, Tc = 2.2695(7), which is in very good agreement 
with the exact value. The FSS ansatz for the (finite lattice) susceptibility maxima 
Xmax is taken as usual as Xmax oc L].^'^. From a (linear) least-square fit, we find that 
7/z/ = 1.75(4) is in perfect agreement with the exact value 7/4. Concerning the specific 
heat we expect in the case of the Onsager exponent a = a logarithmic divergence of 
the form Cmax = a + b\og{Lx). Indeed, the data can be fitted nicely with this ansatz, 
cf. Figure 4. We also tried an unbiased fit using the power-law ansatz Cmax = a + bL'^^'^, 
which gives us a = 0.05(2), verifying the expected value. 
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Figure 4. FSS of the specific- heat maxima Cmax- The logarithmic fit C, 
a + b\og{Lx) and an unbiased fit using the power-law ansatz Cmax = i 
almost indistinguishable on the scale of the figure. 



bLx^'^ are 



For the aspect ratio ( = 0.2, we also sampled the spin-spin correlation functions 
{(JiyCTxy) and {(JxiCTxy) in x- and y-direction, respectively. In the ordered regime where 
nearly all spins are aligned in the same direction {axy = +1, say) we find for {aiyaxy) 
and 1 < y < Ly/2 (where hy = H > 0) almost constant values near unity as one expects. 
The fast decay of the spin-spin correlation function to a slightly smaller value in the 
upper half of the system indicates that the interface is localized near the boundary, cf. 
Figure 5 (a). In the regime with an interface in the bulk along the x-direction we find 
a symmetric shape of {cFiyC^xy) as a function of y which is a clear signal for the phase 
separation, see Figure 5 (b). The opening angle between the plus and minus phases for 
1 < 2/ < Ly/2 and Ly/2 + 1 < y < Ly, respectively, is a measure for the fiuctuations 
of the interface, e.g. a stiff interface shows an acute angle. For temperatures above the 
critical temperature, i.e. in the disordered phase, we observe a similar vanishing of the 
spin-spin correlation functions as in the pure 2D Ising model, see Figure 5 (c). 

At zero temperature, for aspect ratios ( larger than the critical ratio Cs = 0.25 
and strong fields H > J(l -|- 4/Ly), the interface is localized on the boundary, cf. 
Figure 1. Although for > no real solution of Equation (2) exists near Tc, one 
can solve the equation for small temperatures:}: and finds in the case C = 0.5 the phase 
diagram shown in Figure 6. The zero temperature limit is consistent with the energetic 
arguments in Ref. [9], see also Figure 1. The lower line starting at T = 0, H = 1 
can also be detected by means of computer simulations, but as one can see in the 
right plot of Figure 6, for T = 1 the dip between the two maxima in the boundary 
magnetization density mb = {l/Ly)Y.yLiO'i,y vanishes with increasing lattice sizes. 

I Here we are in disagreement with Ref. [9] where no solution was found, because of a mistake in the 
discriminant of Equation (2). 
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Figure 5. Spin-spin correlation function (aiyCTxy) for C = 0.2, Ly — 100, and H = 0.5 
measured in different regimes: (a) ordered regime at T = 1.0, (b) regime with an 
interface in the bulk at T = 2.0, and (c) disordered regime at T = 2.5. 

Therefore, there is no signal for a first-order transition between these two regimes. 
One can argue that there is no phase transition at all, because in the infinite-volume 
limit only the ordered phase survives. The dashed line in the left plot of Figure 6 
starting at T = 0, if = = 2 is not visible in simulations, because this second solution 
of Equation (2) would correspond to the boundary between the bulk magnetized state 
and configurations with an interface propagating inside the bulk which, however, have 
a higher free energy then configurations with an interface localized near the boundary 
and hence are suppressed. Therefore the phase diagram for ( > 0.25 consists only 
of two phases, namely the ordered low-temperature phase and the disordered high- 
temperature phase as in the pure 2D Ising model. We hence conclude that for C > 0.25 
the inhomogeneous boundary magnetic field only leads to finite-size effects. 

Finally, let us come to the special case of ( = (s = 0.25. While the transition line 
disappears for ^ > 1/4 as the solutions of Equation (2) move to the complex plane, for 
(s we still do see different peaks associated with the two phases and, therefore, a finite 
interface tension, cf. Table 2. The numerically estimated transition points for if = 0.5 
and T = 1.5 also contained in Table 2 are again seen to be in good agreement with 
Equation (2). Furthermore, we analysed for T = 1.0 the spin-spin correlation functions 
WiyCxy) and {cTxicr^y) in x- and y-direction, respectively. Slightly below the transition 
from the ordered phase with (m) = mo to the phase with an interface propagating 
inside the bulk and therefore (m) = 0, {aiya^y) shows an asymmetric shape and a fast 
decay near the boundary, indicating that the interface is localized near the boundary, 
cf. Figure 7 (a). With increasing boundary magnetic field we cross the first-order line at 
Hq ~ 0.925, where the interface starts moving into the bulk and, therefore, the profile 
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Figure 6. Left plot: The finite-size phase diagram for a system with C = 0.5. The 
thick hne shows the transition hne between the ordered phase and the region with an 
interface localized on the boundary. The dashed line indicates where a configuration 
with no interface and a configuration with an interface propagating inside the bulk have 
the same free energy (non-physical solution of Equation (2) , because this bulk interface 
configuration has a higher free energy then the one with an interface localized on the 
boundary and hence is suppressed.). The thin vertical line indicates the second-order 
bulk phase transition. Right plot: Probability density of the boundary magnetization 
TOf, for T ™ 1.0 plotted for various lattice sizes ranging from Ly — 16 to Ly — 88 and 
boundary magnetic field values H ^ I where the peaks are of equal height. 

of {(TiyCxy) becomes symmetric in y. Right at the transition hne where we have two 
coexisting phases we find both the fast decay near the boundary as well as the opening 
angle between the two almost symmetric halves, see Figure 7 (b). When the boundary 
magnetic field is increased further, this mixed-phase effect vanishes and the interface 
in the bulk becomes stable. In this case we find a symmetric shape of the spin-spin 
correlation function, cf. Figure 7 (c), similar to C, = 0.2 in Figure 5 (b). 

4. Summary 

Our Monte Carlo data clearly confirm the theoretical considerations of Clusel and 
Fortin [9] and extend their exact results by studying the cases C equal and larger than 
the critical value C,s = 1/4. The observed finite-size scaling behaviour fits nicely with 
their predictions for the infinite system, cf. our results in Table 2. We also find that 
for a large aspect ratio some interesting finite-size effects can be observed, such as, for 
example, a regime in the H-T plane with two states separated by an energy gap which 
vanishes in the infinite-volume limit. Furthermore, we studied the spin-spin correlation 
function for which analytical results are not yet available. Since this observable turned 
out to be quite sensitive to the interface location, it would be a challenging enterprise 
to pursue further analytical considerations in this direction. 
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Figure 7. Spin-spin correlation function (aiyUxy) for C = 0.25, Ly — 80, and 
T = 1 measured for different values of the boundary magnetic field: (a) H — 0.9, 
(b) H = 0.925, and (c) H = 0.975. 
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